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ABSTRACT 

We calculate merger rates of dark matter haloes using the Extended Press- 
Schechter approximation (EPS) for the Spherical Collapse (SC) and the Ellip- 
soidal Collapse (EC) models. 

Merger rates have been calculated for masses in the range lO lo M h _1 to 
lO 14 M h _1 and for redshifts z in the range to 3 and they have been com- 
pared with merger rates that have been proposed by other authors as fits to the 
results of N-body simulations. The detailed comparison presented here shows 
that the agreement between the analytical models and N-body simulations de- 
pends crucially on the mass of the descendant halo. For some range of masses 
and redshifts either SC or EC models approximate satisfactory the results of N- 
body simulations but for other cases both models are less satisfactory or even bad 
approximations. We showed, by studying the parameters of the problem that a 
disagreement -if it appears- does not depend on the values of the parameters 
but on the kind of the particular solution used for the distribution of progenitors 
or on the nature of EPS methods. 

Further studies could help to improve our understanding about the physical pro- 
cesses during the formation of dark matter haloes. 

Subject headings: galaxies: halos - formation -structure; methods: numerical 
-analytical; cosmology: dark matter 
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1. Introduction 

The development of analytical or semi-numerical methods for the problem of structure 
formation in the universe helps to improve our understanding of important physical 



processes. A class of such methods is based on the ideas of Press & Schechter (1974) and 



on their extensions (Extended Press-Schechter Methods EPS, Bond et al. (1991), Lacey 



& Cole (1993)): The linear overdensity <5(x; R) = [p(x; R) — p&(x; i?)]/p fe (x; R) at a given 
point x of an initial snapshot of the Universe fluctuates when the smoothing scale R 
decreases. In the above relation, p(x; R) is the density at point x of the initial Universe 
smoothed by a window function with smoothing scale R. The index b denotes the density 
of the background model of the Universe. This fluctuation is a Markovian process when 
the smoothing is performed using a top-hat window in Fourier space. For any value of 
the smoothing scale R, the overdensity field is assumed to be Gaussian with zero mean 
value. The dispersion of these Gaussians is a decreasing function of the smoothing scale R 
reflecting the large scale homogeneity of the Universe. The mass M contained in a given 
scale R depends on the window function used. For the top-hat window this relation is: 
M = ^7Tp b i R 3 = -R 3 , where p^ and {l m ^ are the values of the mean density and the 

density parameter of the Universe, G is the gravitational constant and Hi is the Hubble's 
constant. The index i indicates that all the above values are calculated at the initial 
snapshot. The dispersion in mass, a 2 , at scale R is a function of mass M and it is usually 
denoted by S, that is S(M) = a 2 [R(M)]. 

In the plane (S, 5) random walks start from the point (S = 0, 5 = 0) and diffuse as S 
increases. Let the line B = Bsc(z) that is a function of redshift z. In the case this line is 
parallel to S'-axis in the (S, 5) plane, then it has a physical meaning as it can be connected 
to the spherical collapse model (SC): It is well known that in an Einstein-de Sitter Universe, 
a spherical overdensity collapses at z if the linear extrapolation of its value up to the present 



exceeds 5 SC ~ 1.686 (see for example Peebles (1980)). All involved quantities (density 
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overdensities, dispersion) are linearly extrapolated to the present and thus the barrier in 
the spherical collapse model is written in the from B(z) = 1.686/ D(z), where D(z) is the 
growth factor derived by the linear theory, normalized to unity at the present epoch. It is 
clear that the line B(z) is an increasing function of z. If a random walk crosses this barrier 
for first time at some value So of S, then the mass element associated with the random 
walk is considered to belong to a halo of mass Mq = S~ 1 (Sq) at the epoch with redshift 
z. However, the distribution of haloes mass Jm, at some epoch z, is connected to the first 
crossing distribution fg, by the random walks, of the barrier that corresponds to epoch z 
with the relation: 

f M (M)dM = f s (S)\^^\dM (1) 

A form of the barrier that results to a mass function that is in better agreement with the 
results of N-body simulations than the spherical model is the one given by the Eq. 

B EC (S,z) = ^B sc (z)[l + (3[S/aB 2 sc (z)n (2) 

In the above Eq. a, (3 and 7 are constants. The above barrier represents an ellipsoidal 



collapse model (EC) (Sheth & Tormen 1999). The barrier depends on the mass 



(S = S(M)) and it is called a moving barrier. The values of the parameters are 
a = 0.707, (3 = 0.485, 7 = 0.615 and are adopted either from the dynamics of ellipsoidal 
collapse or from fits to the results of N-body simulations The spherical collapse model 
results for a = 1 and (3 = 0. 

In a hierarchical scenario of the formation of haloes, the following question is fundamental: 
Given that at some redshift zq a mass element belongs to a halo of mass Mo, what is the 
probability the same mass element at some larger redshift z (z > zq) -that corresponds to 
an earlier time- was part of a halo with mass M with M < M ? This question in terms 
of first crossing distributions and barriers can be written in the following equivalent form: 
Given a random walk passes for the first time from the point (8q, Sq) what is the probability 
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this random walk crosses a barrier B with B > 5q, for the first time between S, S + dS 
with S > S 1 



If we denote the above probability by f(S/8o, So)dS it can be proved, (Zhang & Hui 2006), 
that for an arbitrary barrier, / satisfies the following integral equation: 



f(S/S , S ) = gi(S, 5 , So) + / g 2 (S, S')f(S'/6 , S )dS' 



where: 



gi(S,5 ,S 



92(S,S') 



B(S)-5p 2 dB(S) 
S — 5*o dS 
dB(S) B(S) - B(S' 



dS 



S-S> 



P [B(S),S/5o,So] 
P [B(S),S/B(S'),S' 



and 



P (x,y/x ,y Q ) 



1 



-e ' 2A y 



(3) 

(4) 
(5) 

(6) 



VMKy) 
with Ay = x — Xo and Ay —y — yo- 

In the case of a linear barrier Eq.(3) admits an analytic solution. If B(S) = u + qS, 
where the coefficients u and q could be functions of the redshift z in order to describe the 
dependence on the time, the solution is written: 



f(S/6 ,S 



B(S ) - 5o 



: CXp 



[BjS) - S Q \ 
2(S - So) 



(7) 



'27r(S - So) 3 

Thus, the spherical model which is of the form B = B(z) = u(z) = 1.686/ D(z) leads to the 
solution: 

(Ac) 2 



fsc(S, Zj S , Zq) 



Aoj 



exp 



2AS 



v /2tt(A5)3 

where, AS = S — S , and Aw = cu(z) — co(z ) . 

Unfortunately, no analytical solution exists for the ellipsoidal model. The exact numerical 



solution of Eq.(3) is well approximated by the expression proposed by Sheth & Tormen 



(2002) that is: 



, fQ /Q x 1 \T(S,z/S ,Zq)\ 

f EC {S,z/So,zo) = -j= ^ Ag p /2 ex P 



(Ag) 
IAS 



21 



dS 



(9) 
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where, AB = B(S, z) — B(Sq, z ), and the function T is given by: 

5 \q _ q]n an 

T(S, z/S , z ) = B(S, z) - B(S , z ) + n! 8& B ( S > *)■ ^ 

n=l 

According to the hierarchical clustering any halo is formed by smaller haloes (progenitors). 
A number of progenitors merge at z and form a larger halo of mass M at Zq (zq < z). 
Obviously, the sum of the masses of the progenitors equals to M . Given a halo of mass M Q 
at zq the average number of its progenitors in the mass interval [M, M + dM] present at z 
with z > Zq is : 

^(M/M ,Au)dM= ^f(S,z/S ,z )dM (11) 

Recent comparisons show that the use of EC model improves the agreement between 
the results of EPS methods and those of N-body simulations. For example, [Yahagi et 



al. (2004) showed that the multiplicity function resulting from N-body simulations is far 



from the predictions of spherical model while it shows an excellent agreement with the 



results of the EC model. On the other hand, Lin et al. (2003) compared the distribution 
of formation times of haloes formed in N-body simulations with the formation times of 
haloes formed in terms of the spherical collapse model of the EPS theory. They found 



that N-body simulations give smaller formation times. Hiotelis & Del Popolo (2006) 
showed that using the EC model, formation times are shifted to smaller values than those 
predicted by a spherical collapse model. Additionally, the EC model combined with the 



stable "clustering hypothesis" has been used by (Hiotelis 2006) in order to study density 
profiles of dark matter haloes. Interesting enough,, the resulting density profiles at the 
central regions are closer to the results of observations than are the results of N-body 
simulations. Consequently, the EC model is a significant improvement of the spherical 
model and therefore we are well motivated to study merger-rates of dark matter haloes for 
both the SC and the EC model. This study depends upon the accurate construction of a 
set of progenitors for any halo for a very small "time step" Au. The set of progenitors are 
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created using the method proposed by Neinstein & Dekel (2008) that we describe in Sect. 3. 
In Sect. 2 we define merger rates and we recall fitting formulae resulting from N-body 
simulations. In Sect. 4 our results are presented and discussed. 



2. Definition of merger rates and analytical formulae 

We examine descendant haloes from a sample of Nd haloes with masses in the range 
Md, Md + dMd present at redshift za- For a single descendant halo the procedure is as 
follows: Let M P) i, M Pj 2...M Pt f. be the masses of its k progenitors at redshift z p > Zd- For 
matter of simplicity we assume that the most massive progenitor is M Pj \. We define 
& = M Py i/M P) i for i > 2 and we assume that the descendant halo is formed by the following 
procedure: During the interval dz = z p — Zd every one of the progenitors with i > 2 merge 
with the most massive progenitor i = 1 and form the descendant halo we examine. We 
repeat the above procedure for all haloes in the range Md, Md + dMd found in a volume V 
of the Universe. Then, we find the number denoted by N of all progenitors with £j, i > 2 
in the range (£, £ + d£) and we calculate the ratio N/(VdzdM d d^). We define the merger 
rate B m as follows: 

N 

B ^ M ^^ : ^=vd^m (12) 

Let the number density of haloes with masses in the range Md, Md + dMd at Zd be 
n{M d ,Zd) = Nd vdMa d ^ ■ T ne ra ti° B m /n = N/(N d dzd^) measures the mean number of 
mergers per halo, per unit redshift, for descendant haloes in the range Md, Md + dMd with 
progenitor mass ratio £. 



Fakhouri & Ma (2008) analyzed the results of the Millennium simulation of Springel et 



al. (2005). The fitting formula proposed by the above authors is separable in the three 



variables, mass Md, progenitor ratio £ and redshift z: 
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B(M d ,£,z p : z d ) 
n(M d , z) 



A ■ F{M d )G{i)H{z) 



(13) 



with 

F(M d ) 



M 



«i 



G(0 = r 2 exp 
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H(z) = (^) 4 and the values of the 

-2.01,0a 



parameters are M = 1.2 x 1O 12 M ,A = 0.0289, £ = 0.098, a x = 0.083, a 2 
0.409, a 4 = 0.371. 



Lacey & Cole ( 1993 ) showed that in the spherical model the transition rate is given by: 

dcr(M) 



r ( M — ► M d /z d )dM d = (2/tt) 1/2 



d5 c {z) 
dz 



\ o\M d y 


-3/2 




(T 2 (M) _ 


exp 





1 



a 2 (M d 



1 



v cr 2 (M d ) a 2 (M)J_ 



dM 
dM d 



M=M d 



(14) 



This provides the fraction of the mass belonging to haloes of mass M that merge 
instantaneously to form haloes of mass in the range M d , M d + dM d at z d . The product 
r ■ f sc (M, z d )dM, where f sc (M,z) is the unconditional first crossing distribution for the 
spherical model, gives the above fraction of mass as a fraction of the total mass of the 
Universe and successively multiplying by (pb/M) ■ V the number of those haloes is found. 
Then, dividing by (pb/M d ) ■ V ■ f sc (M d , z d )dM d (that equals to the number of the descendant 
haloes) we find: 



N 



N,dz 



d5 c (z) 
dz 



Ma 



M a 2 (M) 



do-(M) 
dM 



o\M d 



M=M d 
3/2 

dM 



(15) 



CT 2 (M) _ 

Assuming a strictly binary merger history i.e. every halo has two progenitors, and denoting 
by £ the mass ratio of the small progenitor to the large one (£ = (M d — M)/M), using 
dM = |f^d£ and substituting in (15) we have the final expression for the binary spherical 
case, that is: 



B 



in 
n 



N 
Nddzdi 



dd c (z) 
dz 



M 
a 2 {M) 



da(M) 
dM 



J M=M d 
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o\M d ) 
a 2 (M) 



-3/2 



(16) 



3. Construction of the set of progenitors 



The construction of progenitors of a halo can be based either or Eqs (8) and (9) or else 
on Eq. (11). For the first case a procedure is as follows: A halo of mass M at redshift z$ 
is considered. A new redshift z is chosen. Then, a value AS is chosen from the desired 
distribution given by Eq.(8) or (9). The mass M p of a progenitor is found by solving for 
M p the equation AS = S(M P ) — S(M ). If the mass left to be resolved M — M p is large 
enough (larger than a threshold), the above procedure is repeated so a distribution of the 
progenitors of the halo is created at z. If the mass left to be resolved -that equals to Mq 
minus the sum of the masses of its progenitors- is less than the threshold, then we proceed 
to the next time step , and re-analyze using the same procedure. 



A complete description of the above numerical method is given in Hiotelis fc Del Popolo 



(2006). The algorithm - known as N-branch merger-tree- is based on the pioneer works of 



Lacey & Cole (1993), Somerville & Kollat (1999) and van den Bosch (2002) 



We have to note that the construction of a set of progenitors for an initial set of haloes 
after a "time step" Au is a problem that has not a unique solution. Consequently, it is 
interesting to compare different solutions with the results of N-body simulations in order 
to find those which show a better agreement. We note that any of the above proposed 
algorithms has a number of drawbacks. The algorithm to be used has to be suitable for the 
particular problem. If for example the algorithm assumes an initial set of descendant haloes 
of the same mass, it cannot be used for more than one time steps since the set of progenitors 
predicted at the first time step does not consist of haloes of the same mass. Since our 
purpose is the derivation of merger rates, we used the method proposed by Neinstein fc 



Dekel (2008) that is suitable for the calculation of a set of progenitors for descendant haloes 
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of the same mass for a single time step. A description is given below: 
We assume a set of N d haloes of the same mass M at z = z d . We use the variables 
Mi, M 2 , M 3 .. to denote the masses of their progenitors at redshift z p , after a time step 
Alj = ui(zd) — uj(z p ). We assume that M i > M 2 > M 3 , ... and we denote by Pi(M) the 
probability that the i th progenitor has mass M. We also assume that the value of M 1; that 
is the mass of the most massive progenitor of a halo, defines with a unique way the masses 
of all its rest progenitors. Additionally, Pj/i(Mj/Mi) is the constrained probability that the 
jth p ro g en itor of a halo equals M given that its most massive progenitor is M x . Obviously 
the following Eqs. hold: 

P t (M) = J P t/ i(M/Mi)Pi(Mi)dMi (17) 
P tot (M) = J2^( M ) ( 18 ) 

i 

P(M u M 2 ,...) = if ^M,>M (19) 

i 

These are the key equations for the construction of the set of progenitors. We use the 
following three steps: 

1st step: The distribution of the most massive progenitors. 
We define P tot (M) using Eq.(ll), that is: 

P tot (M)dM = Sm/Mq, Auj)dM (20) 
aM 

The value of the integral P to t(M)dM depends on both Mum and Aco. For M min — > 

it declines due to the presence of the large number of very small progenitors. The value of 
the integral increases for increasing Au. Thus, for reasonable choice of M min the values of 
the above integral is larger than unity. Then, the distribution of Mi can be found by the 
following procedure: First, we solve the Eq. 

/ P t ot(M)dM = 1 (21) 

J x*Mq 
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with respect to x*. The resulting values of M* = x*M are larger than M min . Then, we pick 
Mi from the distribution: 

f Ptoi(Mi), Mi > M* 

Pim = { (22) 

I 0, otherwise 

This is done by the following procedure: A random number r is chosen in the interval [0, 1] 
and the equation P tot (M)dM = r is solved for Mi. The resulting values of Mi have the 
above described distribution. 

If Mi e ft = Mo — Mi > M m i n we proceed with the second progenitor. Otherwise, the halo 
has just one progenitor and we proceed with the next halo. 
2st step: The distribution of M 2 . 

Let fi(Mi) be the mass of the i th progenitor given that the mass of the most massive 
progenitor equals to Mi. We assume that 

P l/1 (M t /M 1 ) = 5[M t - /-(Mi)] (23) 

where 5 is a delta function and /, a monotonically decreasing function of Mi. 
We consider the differential equation: 

d/ i (M 1 ) Pi(Mi) 



dMi " PMMi)} 
Using (23) and (24) the right hand side of (17) is written: 

5[Mi - /.(MOJP^MOdMi = 



(24) 



6[Mi - MMjPiiMM^dMMi) = / 6[Mi - UiM^P^M^f^) = 

J — oo 

PMMi)} = Pi(Mi) (25) 

and thus the solution of the differential Eq. (24) satisfies Eq. (17). Thus, the mass of the 
second progenitor can be found by integrating numerically (24) for i = 2: 

d/ 2 (Mi) Pi(Mi) 



dMi P 2 [/ 2 (Mi)] 



(26) 
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' ' ' ' I ' ' ' ' I ' ' ' ' I ' ' ' ' I ' ' ' ' I ' ' ' ' I ' ' ' ' I ' ' ' ' I ' ' ' ' I ' ' ' 
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 

M/M n 



Fig. 1. — The distribution of the first and the second progenitors M x and M 2 , respectively 
for both the SC and EC models. Filled gradients show the distribution of Mi and empty 
gradients show the distribution of M 2 for the spherical model for M = 100, Acu = 0.1 and 
Zd = 0. Squares show the same distributions for the ellipsoidal model. Solid and dashed 
lines are the predictions of Eq.(ll) with / given by (8) and (9), respectively. 
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Fig. 2. — The function f 2 (M 1 ), the solution of the differential equation (26), equals the mass 
M 2 of the second progenitor and is plotted for M = 100 and Auj = 0.1. For the same 
mass Mi of the most massive progenitor M 2 is smaller for EC model than for SC model. 
Additionally, the ellipsoidal model extends to lower values than the spherical model. The 
lowest values of Mi shown, are about 0.36M for the ellipsoidal model and 0.41M for the 
spherical model, respectively. 
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The function P 2 involved is unknown. So a trial function ^2(3/) = Ptot{y) ~ Pi(y) is used 
and the Eq. 

% = G{*,y) (27) 

where G(x,y) = — P\(x)/P2(y) is solved numerically for y using a classical 4 th order 
Runge-Kutta with initial conditions Xj n = M*,y in = M 2j o = M* (called solution I in 



Neinstein fc Dekel (2008)). We used a step Ax = [Mi - M*]/N S} where N s defines the 



number of steps. We used various values of N s from 100 to 10000 and we found that the 
results are essentially the same. 

In the case the solution of the above differential equation is M 2 < M m j n then we enforce 
M2 = Mq — Mi. Finally, the resulting values of M2 are used for the numerical construction 
of P 2 . 

In our calculations, we used a flat model for the Universe with present day density 
parameters fl m ,o = 0.3 and Qa,o = A/3-f/g 2 — 0-7- A is the cosmological constant and Hq is 
the present day value of Hubble's constant. We used the value H = 100hKms" 1 Mpc _1 and 
a system of units with m unit = 1O 12 M /i _1 , r unit = l/i _1 Mpc and a gravitational constant 
G — 1. At this system of units, H /H unit = 1.5276. 



As regards the power spectrum, we used the AC DM form proposed by Smith et al. (1998). 
The power spectrum is smoothed using the top-hat window function and is normalized for 
a 8 = a(R = S/i^Mpc) = 0.9. 

We used a number N res = 10 5 haloes of the same mass Mq at zq = Zd and we found their 
progenitors at z p that is after a "time-step" Aw = u(z p ) —u(z d ). We studied three values of 
Zd that are z d = 0,1 and 3 respectively. We examined values of M in the range 0.01 to 100 
in our system of units. These values correspond to masses in the range Mq = lO lo M /i _1 to 
M = 1O 14 M /i _1 . We studied three values of Au namely 0.1,0.05 and 0.025. We also used 
M min = 10" 3 M for Aoo = 0.1 and M min = 5 ■ 10~ 4 M for Aoo = 0.05 and 0.025. 
Fig.l compares the distributions of progenitors for Mq = 100, Au = 0.1 and z d = 0.0 with 
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the analytical ones given by Eq. (11) for both the spherical and the ellipsoidal models. Up 
to this step every halo has at most two progenitors. It is clear that the agreement is very 
satisfactory. 

Fig.2 shows the solution /2(Mi) of Eq. (26). It presents M 2 = /2(Mi) as a function of M 1 
both normalized to M . It corresponds to z d = 3.0 and z p = 3.078, that is to Auj = 0.1. It 
is shown that the distribution of most massive progenitors extends to smaller values in the 
ellipsoidal model. The reflection of this different behavior to merger rates will be studied in 
Sect.4 

The satisfactory agreement between the distributions of progenitors predicted by the 
method studied and by Eq.(ll) holds also for various values of the descendant halo and 
various redshifts. This is shown in Fig.3 where the distribution of progenitors for both SC 
and EC models for z d — 0, 1 and 3 are presented. The value for the "time-step" is Au = 0.1. 
The corresponding values of 0.1097, 1.083 and 3.078, respectively. 

However, if we focus on small values of M/M we see that the distribution of progenitors 
there differs significantly from the theoretical one. Such an example is given in Fig. 4 where 
the thin solid line is the theoretical distribution and the thick solid line is the distribution 
that results after the above two steps. (Dashed line is the final distribution after the 
completeness of 3 d step and it will be discussed later.) This disagreement shows clearly 
that the number of small progenitors is underestimated when every halo is analyzed to two 
progenitors and the need of more progenitors is clear. Although the disagreement appears 
only for small values M/M is important for the calculation of merger rates as it will be 
shown below. 

The above two steps are completed for the whole sample of descendant haloes. Thus, after 
the completion of the second step, the distribution P 2 is found numerically and is expressed 
by a polynomial in order to be used in the 3 d step below. 
3 d step: The distribution of M h i > 2. 
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M/M M/M 



Fig. 3. — The distribution of progenitors cW/dM versus their mass M normalized to the 
mass M of the descendant halo. The first column corresponds to the spherical model for 
Zd — 0, Zd — 1 and = 3 (from top to bottom) and the second column to the ellipsoidal 
model. A value Au = 0.1 is used. Dashed lines are the predictions of the method studied in 
this paper while solid lines are the predictions of Eq. (11) 
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Fig. 4. — The distribution of progenitors dN/dM for small values of M/M for M = 
100, Zd = 0, Auj = 0.1 and for the SC model. Thin solid line is the prediction of Eq. (11) 
and the thick solid line is the prediction of the two first steps of the method studied, that is 
without progenitors Mj with % > 2. Dashed line is the final distribution after the third step, 
that is after the prediction of the full set of progenitors. 
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Obviously, a halo has a progenitor % if the mass left to be analyzed M ie f tti = M — X^n=i _1 -^n 
is Mi e f t>i > M min . We found the distribution of the rest progenitors using the following: 
First, we found the solution R of the equation p2(y) = Po(y)- Obviously for y < R we have 
P 2 (y) < Po(y). Then, we define: 

, * j Pq(x)-P 2 (x), when M min < x < M high;i 
Pi{x) = t (28) 

t 0, otherwise 

where 

M highti = min{Mj e/t)i , M^} for i>3 and M highj3 = R-M (29) 
Finally, we solve Eq. (24) for f^M^. 

Dashed line in Fig. 4 is the distribution of progenitors after the completeness of the third 
step. It is clear that this distribution is much closer to the theoretical one given by the thin 
solid line than the distribution -that is described by the thick solid line- that results using 
only the first two progenitors Mi and M 2 . 



4. Results 

We have already mentioned that distributing progenitors according to Eq. (11) is a 
problem that has not a unique solution. Additionally, the calculation of merger rates of 
dark matter haloes using analytical methods involves a large number of parameters. These 
are: the background cosmology, the model of collapse used (SC or EC), the mass of the 
descendant haloes M , the redshift Zd and the "time step" Au. 

The background cosmology used has been described in the previous section. The 
distribution of progenitors is done according to the method analyzed through this paper. 
So the parameters that were studied are: the model of collapse, the mass of the descendant 
haloes M , the redshift z d and the time step Aui. 
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We give a first result in Fig. 5. Solid line shows the predictions of the formula given by 
Eq.(13), proposed to fit the results of N-body simulation. Thick long dashes show the 
predictions of the spherical binary model given by Eq (16). It is shown that the spherical 
binary model overestimates merge rate for large values of £ while it underestimates the 
merger rate for values of £ smaller than 1CT 2 . The predictions of the method studied are 
shown by squares and thick small dashes. Squares show the results after the first two steps 
described in Sect. 3, that is after the distribution of the two first progenitors Mi and M 2 
only while thick small dashes show the prediction for the whole set of progenitors. The 
third step in the procedure described in Sect. 3 adds progenitors that have small masses. 
This increases the number of progenitors with small £ and rises the curve of the merger 
rate. This result agrees better with the predictions of N-body simulations. The results 
correspond to z d = 0. We used Au = 0.025 that results to z p = 0.02804. 

The accurate calculation of the merger rates requires that Au — > 0. However, we examined 
different values of Au and we verified that the results do not depend crucially on this 
parameter. We used three values of Au; namely 0.025,0.05 and 0.1. Differences in merger 
rates due to the different values of Au are negligible. As an example we present Fig. 6. 
It refers to the SC model for Zd — and for a descendant halo with mass M = 100, for 
Auj = 0.1 and Alu = 0.025 (solid line and dashed line respectively). The corresponding 
values of z v are 0.1097 and 0.02804. Thus Az is about four times smaller in the second case. 
It is clear that only negligible differences are present. 

In Figs 7 and 8 we present results for different masses, redshifts and time-steps. In all 
snapshots dashed lines are the predictions of the SC model and solid lines show the results 
of the EC model. Squares are the predictions of the N-body fitting formula formula given 
by Eq.(13). 

From the results presented in the first row of Fig. 7 is clear that the EC model results 
to merger rates that are not in agreement with the results of N-body simulations, for a 
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Fig. 5. — Merger rate for M = 100 for z d = and z p = 0.02804 (that corresponds to 
Auj = 0.025). Solid line corresponds to the formula proposed to fit the results of N-body 
simulation that is given by Eq.(13). Thick long dashes show the predictions of the spherical 
binary model given by Eq.(16). Squares are the predictions of the method studied in this 
paper for the SC, using only two progenitors Mi and M 2 . Thick small dashes show the 
prediction of the above method for the whole set of progenitors. 
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descendant halo of small mass M = 0.01. Instead the results of SC seem to be satisfactory. 
For larger masses the agreement between EC model and N-body results becomes better. 
For large haloes, M = 100, EC model approximates N-body simulations better than SC 
model. All the results of Fig.7 have been calculated for Aoo = 0.1. Fist column shows 
results for z d = while the second one for z d = 1. 

All curves in Fig.8 have been predicted for Au = 0.025. Three rows correspond to 
M = 0.1, M = 1 and M = 100, respectively. As in Fig.7, different lines represent different 
models. Thick dashes show the results of SC model, solid line the results of EC model and 
squares the results of the fitting formula give by Eq. (13). 

A more detailed comparison between the results of EPS and N-body simulations is given 
in Fig.9. We calculated the relative difference (Reps — Rnb)/Rnb where Reps and R NB 
are the merger rates predicted by the EPS and by N-body results, respectively. The results 
presented in this Fig. can be summarized as follows: 

For low redshifts (z = to z = 1), merger rates of haloes with descendant mass in the 
range lO lo M h _1 to lO 12 M h _1 derived by the SC model fit very satisfactory the results 
of N-body simulations. For example, for z d = and Alu = 1 the difference is less than 15 
percent, except for some very small values of x. Instead, for the same range of masses and 
redshifts, merger rates derived by the EC are significantly lower than those predicted by 
N-body simulations. 

For the above range of redshifts and for haloes of mass 10 14 M o h~ 1 the fits by EC model 
are very satisfactory (in general the relative difference is smaller than 20 percent) while the 
results of SC are significantly higher than those of N-body simulations. 
For a higher redshift (z = 3) both SC and EC model overestimate the merger rate of 
large haloes. Merger rates of smaller haloes are overestimated by the SC model and 
underestimated by the EC model. The above conclusion seems not to depend, at least 
significantly, on the values of redshift z d and time-step Au. 
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We have to note here that both N-body simulations and analytical methods have problems 
in describing very accurately some physical properties of dark matter haloes. This is due 
to either technical difficulties or to the fact that some physical mechanisms are not taken 



into account. For example, in a recent paper Fakhouri & Ma (2010) use the results of the 



Millenium - II simulation, (Boylan et al. 2009), to derive a formula of the same form 
of that given in Eq.(13). Millenium - II simulation has better resolution than Millenium, 



(Springel et al. 2005), simulation. Due to the better resolution, the best fitting values of 
the parameters in Eq.(13) are changed. For example, the value of a\, that is the exponent 
of the mass of the halo, from 0.083 becomes now 0.133. Obviously the dependence of mass 
remains weak but such a change in the value of at results, for a halo with Md = 100, to a 
new merger rate that is 26% larger. This percentage is too large since it can change the 
whole picture, at least for large haloes, resulting from our comparison. Additionally, it is 
interesting to notice Fig. Al in the appendix of the above paper. It describes merger rates 
given by five different algorithms. These algorithms are used to analyze the results of the 
same simulation and to study fragmentation effects in FOF (friends of friends) merger trees. 
From this Fig. it is clear that differences due the use of different algorithms may be larger 
than the differences between analytical methods and N-body simulations derived by our 
study and shown in our Fig. 9. 

From the above discussion it is clear that the results of N-body simulations are very 
sensitive not only to the resolution but also to the halo finding algorithm. This sensitivity 



can lead to completely different results. The following example is very characteristic: Bet et 



al. (2007) studied, among other things, the value of the spin parameter as a function of the 



mass of dark matter haloes. They found that the FOF algorithm results to a spin parameter 
that is an increasing function of mass while a more advanced halo finding algorithm, that 
they been proposed, results to a spin parameter that is a decreasing function of mass! On 
the other hand, N-body simulations have the ability to deal with complex physical process. 
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For example the destruction of dark matter haloes as well as the the role of the environment 
are factors that are not taken into account in most of the analytical methods. This is an 
additional reason for the presence of differences between the results. 
Summarizing our results we could say that: SC approximates better the merger rates of 
small haloes while EC the merger rates of heavy haloes. This is obviously an interesting 
information, but since it has been resulted from a specific solution for the problem of 
the distribution of progenitors, a further study of different solutions is required. The 
finding of a solution that approximates satisfactory merger rates from N-body simulations, 
independently on the redshift and mass should be an important achievement. Such a trial 
requires future comparisons and obviously improvements on both kind of methods. 
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Fig. 6. — The role of the 'time-step' in the estimation of merger rates: For the EC model, 
M = 100 and z d = solid and dashed lines correspond to Aoo = 0.1 and Acu = 0.025, 
respectively. It is clear that differences are negligible. 
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Fig. 7. — Merger rates for various models and various values of M and z&. The two snapshots 
of the first row show merger rates for M = 0.01 at — and Zd = 1, respectively. Squares 
are the predictions of N-body given by Eq. (13). Thick dashes show the results of SC model 
and solid lines the results of EC model by the method used in this paper for Acu = 0.1. 
Snapshots of the second row correspond to M = 1 and those of the third row to M = 100. 
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Fig. 8. — Merger rates for various models and various values of M and z^. The two snapshots 
of the first row show merger rates for M = 0.1 at Zd — and Zd = 3, respectively. Squares 
are the predictions of N-body given by Eq. (13). Thick dashes show the results of SC model 
and solid lines the results of EC model by the method used in this paper for Alu = 0.025. 
Snapshots of the second row correspond to M = 1 and those of the third row to M = 100. 



-27- 



i 

0.8 
0.6 
0.4 
0.2 



-0.2 
-0.4 




N 



Solid line: M =0.01 
0.6 FDashed line: M =l 
_ 8 |Dashdotdot: M A =100 



i/v V\aa_.. 



Spherical Model 
z d =0, A(O=0.1 

-X 1 



Ellipsoidal Model 
z d =0, Aco = 0.1 



'• V \./v, ,\ 



0.8 
0.6 
Z 0.4 
^ 0.2 



Spherical Model 
z=0,Aco =0.025 



Pi 



o 




Ellipsoidal Model 
z d =0, Aco=0.025 



Y*J\hJ \ 



N -./\ 



A 



r - 2 
& -°- 4 

-0.6 
-0.8 
-1 



Solid line: M =0.1 
Dashed line: M =l 
Dashdotdot:M n =100 




0.8 i 



0.4 
0.2 


-0.2 
-0.4 
-0.6 
-0.8 
-1 




Solid line: M =0.1 
Dashed line: M =l 
Dashdotdot: M„=100 



Spherical Model 
z d =3, Aco=0.025 



Ellipsoidal Model 
z d =3, Aco=0.025 



0.25 



0.5 



0.75 



0.25 



0.5 



0.75 



Fig. 9. — Detailed comparisons between the EPS and N-body results. The relative difference 
(Reps — Rnb)/Rnb, where Reps and Rnb are the merger rates predicted by the EPS and 
by N-body results respectively, is plotted as a function of £. SC model gives merger rates 
that are in good agreement with N-body results for small haloes (in the range 0.01 — 1) while 
EC model approximates better the merger rates of heavy haloes (M = 100). 
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